function VARbs = proxy_svar_bs(VAR,nboot,clevel)

% Preliminaries
data_Y = VAR.data.Y;
data_Z = VAR.data.Z;
 
[T, n_y] = size(data_Y);
[T_z, k] = size(data_Z);  
 
p = VAR.p; 
T_y  = T - p;  
[overlap,overlap_idx] = intersect(VAR.data.time_y,VAR.data.time_z,'stable'); 
ty_zstart = overlap_idx(1);
ty_zend = overlap_idx(end);


% Wild Bootstrap
  jj=1;
    % removes a constant
    res = detrend(VAR.res,'constant'); 
     while jj<nboot+1
       %T x 1 -1 or 1 random numbers with prob 0.5-0.5
       rr = 1-2*(rand(VAR.T,1)>0.5);
       %T x n randomly choose the sign of the time T shocks (all)
       resb = (res.*(rr*ones(1,VAR.n)))';
       varsb = zeros(VAR.p+VAR.T,VAR.n); 
       
       % Initial values
       varsb(1:VAR.p,:)=VAR.data.Y(1:VAR.p,:);
        for j=VAR.p+1:VAR.p+VAR.T
            lvars = (varsb(j-1:-1:j-VAR.p,:))';
            varsb(j,:) = lvars(:)'*VAR.beta(1:VAR.p*VAR.n,:)+VAR.beta(VAR.p*VAR.n+1:end,:)+resb(:,j-VAR.p)';
        end 
              
        VARBS = VAR;
        VARBS.data.Y = varsb; 
        VARBS.data.Z = [VAR.data.Z.*(rr(ty_zstart-p:ty_zend-p,1)*ones(1,size(VAR.data.Z,2)))];
        VARBS = proxy_svar(VARBS);

        VARbs.beta(:,:,jj) = VARBS.beta;

        for j=1:VAR.k
            irs = VARBS.irs(:,:,j); 
            VARbs.irs(:,jj,j) = irs(:);
        end  
      jj=jj+1;   

     end
     
% bias correction    
bias = squeeze(mean(VARbs.beta,3)) - VAR.beta;
VARbs.beta = VAR.beta + bias;

% stability check
bb = zeros(VAR.n,VAR.n*VAR.p);
bb(1:VAR.n,1:VAR.n) = VAR.beta(1:VAR.n,1:VAR.n);
for jjj=2:VAR.p
    bb(1:VAR.n,1+(jjj-1)*VAR.n:jjj*VAR.n) = VAR.beta(1+(jjj-1)*VAR.n:jjj*VAR.n,1:VAR.n);
end
d = stability(bb,VAR.p); 
if d >= 1; error('--> WARNING: the estimated VAR model is not stable!\n'); end

% Impulse Responses 
irs(VAR.p+1,:) = VAR.b1;
for jj=2:VAR.IRF_hor
    lvars = (irs(VAR.p+jj-1:-1:jj,:))';
    irs(VAR.p+jj,:) = lvars(:)'*VARbs.beta(1:VAR.p*VAR.n,:);     
end
VARbs.irs = irs(VAR.p+1:end,:); 

% Second (!) Wild Bootstrap
  jj=1;
    % removes a constant
    res = detrend(VAR.res,'constant'); 
     while jj<nboot+1
       %T x 1 -1 or 1 random numbers with prob 0.5-0.5
       rr = 1-2*(rand(VAR.T,1)>0.5);
       %T x n randomly choose the sign of the time T shocks (all)
       resb = (res.*(rr*ones(1,VAR.n)))';
       varsb = zeros(VAR.p+VAR.T,VAR.n); 
       
       % Initial values
       varsb(1:VAR.p,:)=VAR.data.Y(1:VAR.p,:);
        for j=VAR.p+1:VAR.p+VAR.T
            lvars = (varsb(j-1:-1:j-VAR.p,:))';
            varsb(j,:) = lvars(:)'*VARbs.beta(1:VAR.p*VAR.n,:)+VARbs.beta(VAR.p*VAR.n+1:end,:)+resb(:,j-VAR.p)';
        end 
              
        VARBS = VAR;
        VARBS.data.Y = varsb; 
        VARBS.data.Z = [VAR.data.Z.*(rr(ty_zstart-p:ty_zend-p,1)*ones(1,size(VAR.data.Z,2)))];

        try
            VARBS = proxy_svar(VARBS, bias);
        catch
            continue
        end

        
        for j=1:VAR.k
            irs = VARBS.irs(:,:,j); 
            VARbs.irs_bs(:,jj,j) = irs(:);
        end  
      jj=jj+1;   

     end

for j=1:VAR.k  
    se = reshape(std(VARbs.irs_bs(:,:,j),[],2),VAR.IRF_hor,size(VAR.irs,2));
    VARbs.irsL(:,:,j)= VARbs.irs + norminv((1-clevel/100)/2) * se;
    VARbs.irsH(:,:,j)= VARbs.irs - norminv((1-clevel/100)/2) * se;
%     VARbs.irsL(:,:,j)=reshape(quantile(VARbs.irs_bs(:,:,j)',(1-clevel/100)/2),VAR.IRF_hor,size(VARBS.irs,2));
%     VARbs.irsH(:,:,j)=reshape(quantile(VARbs.irs_bs(:,:,j)',1-(1-clevel/100)/2),VAR.IRF_hor,size(VARBS.irs,2)); 
end

end